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The effect of on-site electron-electron repulsion U in a band insulator is explored for a bilayer 
Hubbard Hamiltonian with opposite sign hopping in the two sheets. The ground state phase di- 
agram is determined at half-filling in the plane of U and the interplanar hybridization V through 
a computation of the antiferromagnetic (AF) structure factor, local moments, single particle and 
spin wave spectra, and spin correlations. Unlike the case of the ionic Hubbard model which has a 
closely related noninteracting dispersion relation, no evidence is found for a metallic phase inter- 
vening between the Mott and band insulators. Instead, upon increase of U at large V, the behavior 
of the local moments and of single-particle spectra give quantitative evidence of a crossover to a 
Kondo insulator state preceeding the onset of magnetic order. Our conclusions generalize those of 
single-site dynamical mean field theory, and show that including interlayer correlations results in an 
increase of the single particle gap with U. 

PACS numbers: 71.10.Hf, 71.27.+a, 02.70.Uu 
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Introduction: Whether interactions might drive metallic 
behavior in two dimensional disordered systems, where 
disorder just marginally succeeds in localizing all the 
cigenstatcs, is a question that has been the subject of 
considerable experimental and theoretical scrutiny [H-Q. 
It is natural to ask the same question concerning band 
insulators which likewise have vanishing dc conductivity 
in the absence of interactions. In the case of band in- 
sulators, carrier density plays an especially central role, 
since the band must be precisely filled. This lends an ad- 
ditional complexity to the issue, since interactions might 
also give rise to Mott insulating and magnetic behavior. 

The possible connection between disordered interact- 
ing systems, and correlated band insulators is made more 
concrete by considering the Anderson model, where ran- 
dom site energies couple to the local density, and the 
'ionic' Hubbard model (IHM)[3| which has a superlattice 
potential where the site energy has a regular structure, 
taking two distinct values on the sublattices of a bipar- 
tite lattice. On the one hand, it is plausible that the 
same physical effects that could cause a metallic tran- 
sition for random site energies, the reduction of charge 
inhomogencity and resulting derealization of the elec- 
tronic wave functions by interparticle repulsion, would 
also be operative in the patterned case. On the other 
hand, momentum is still a good quantum number in the 
presence of a regular array of site energies, suggesting 
possible differences between the effect of U in the two 
situations. 

The approximations made in the most simple, single 
site, Dynamical Mean Field Theory (DMFT) approach 
to the treatment of electron-electron interaction empha- 
size some of the possible nuances in attempting to elu- 
cidate the physics of correlated band insulators. Single 
site DMFT can capture the band insulator (and how it 



differs from an Anderson insulator) by incorporating a 
density of states with N(Ep) = 0. However, it also min- 
imizes the role of momentum, and hence blurs some of 
the distinction between band and Anderson insulators. 

DMFT has, in fact, been used to explore whether cor- 
relations can drive a band insulator metallic. Garg et 
al. found Q that for the IHM, treated within single site 
DMFT, the band gap becomes zero at a critical U c i, with 
a Mott gap re-emerging at a larger U C 2- In between, the 
system is metallic. A subsequent cluster DMFT study 
of Kancharla et aL0] which incorporated antiferromag- 
netic correlations, found a phase diagram with somewhat 
different topology, but still exhibiting an intermediate 
region which was suggested to have bond ordered wave 
character. 

Model: In this paper, we shall consider a bilayer Hubbard 
Hamiltonian: 

^ = _ *I ( C L,<T C k,;,(T + C k,/,cr C j,/,cr) 

1„ 1, 



(1) 



which provides a specific realization of the effect 
of electronic correlation in band insulators. In 
Eq. [TJ Cj i CT (Cj l a ) are creation(destruction) operators for 
fermions of spin a on site j of layer / = 1,2. The in- 
tralayer hoppings are t, between near- neighbor sites j,k 
of a two dimensional square lattice, and the interlayer 
hopping is V. Correlation is introduced in the model 
through the on-site repulsion U. We have included a 
chemical potential \ii for generality. However, here we 
focus on the half- filled case \ii = 0. 
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If the in-plane hoppings[l(J are chosen with oppo- 
site sign, ii = — t?. = t, Eq. [I] is a band insulator at 

U = with i? q = zty/eq + F 2 . This bears a strong 
resemblance to the dispersion relation of the 2D IHM 
whose superlattice potential A^j(— l) J nj similarly has 
the effect of opening a band gap at q = (tt,tt), altering 

the A = dispersion relation e q to Eq = ±^/e 2 + A 2 . 

Real space Quantum Monte Carlo (QMC)0 has sup- 
plemented DMFT studies of the IHM@ by allowing for 
magnetic order and concluded that U could cause the 
appearance of a metallic phase. However, despite the 
similarity in Eq between the IHM and Eq. [I] the physics 
of the two models is fundamentally different: the bilayer 
has twice as many allowed q points and uniform average 
density, (n^i M ) = i. This is in contrast to the staggered 
charge density wave pattern in the presence of the super- 
lattice potential in the IHM. This difference, combined 
with the local character of the interaction, is at the ori- 
gin of the contrast in the ground state properties which 
we will present. 

As with the IHM, the Hamiltonian in Eq. [1] has been 
previously studied within the DMFT formalism Q , with 
a model hybridizing two bands with identical, semi- 
elliptical density of states. DMFT finds a scenario re- 
markably similar to the one observed for the metal- 
insulator case: a first order transition between band and 
Mott insulator characterized by a discontinuous change 
in the double occupancy. Remarkably, upon increas- 
ing the interaction, DMFT also predicts that the single- 
particle gap should monotonically shrink in stark con- 
trast with the behavior in the Mott phase where the gap 
grows monotonically with U. 

We will explore these issues using determinant Quan- 
tum Monte Carlo (DQMC) [H. This method allows an 



exact calculation [l^ of the properties associated to the 
Hamiltonian in Eq. [1] on lattices of finite spatial extent. 
Here we will show results for systems consisting of two 
sheets of up to N = 14 x 14 sites. 

Magnetic transition: We start our discussion of the 
phases of Eq. [T]by looking for possible antiferromagnetic 
long-range order (LRO). The most direct signature is the 
thermodynamic extrapolation of the in-layer structure 
factor (which has the same value on the two layers), 
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a 3 = C it C iV 



converging to a non-zero value as N — > oo. We will focus 
on antiferromagnetism, 3 = S(ir,ir), which is the ex- 
pected dominant magnetic instability at half-filling. Be- 
cause of the continuous spin symmetry and the fact that 
we are in two dimensions, we expect LRO only at T = 0. 

The finite size scaling analysis necessary to locate the 
critical value of U for onset of magnetic order is presented 



1 


1 1 


- AF Mott Insulator 


t i - 




/o 1 


% 


/ 1 




s °o j~ 












Band Insulator , 


l 


l 1 



O dm/dU 
A Saf 
■ ■ ■ MFT 



O U = 0.0 

O U = 2.0 

A U = 3.7 

□ U = 3.9 
U = 4.0 

< U = 4.2 

> U = 4.4 



0.8 
V 

> > 



IS 
a 



V = 0.5. |i- HI. St - 1/8 
I 



0.06 - 
0.04 - 
0.02 





FIG. 1: Finite-size scaling of the AF structure factor at 
V = 0.5. Symbols are DQMC data for S af . Lines are fits 
to third-order polynomials in the inverse linear lattice size 
1/yfN. The inset shows the transition line computed with 
DQMC (dashed) and MFT (dotted). The vertical line corre- 
sponds to the critical V , predicted by studies on the Heisen- 
berg model, above which no magnetic long order is possible. 
Circles correspond to maxima in dm/dU as a function of U 
at constant V (see Fig. 



in Fig. [T]for V — 0.5. At weak coupling U < 4.0 there is 
no LRO. As the on-site repulsion increases, LRO sets in 
around U « 4.2. Our value of U c is significantlysmaller 
than the DMFT estimate of approximately 5.5 [8( which 
was however computed for a transition to a Mott insu- 
lating state without long range magnetic order. 

We repeated the finite size scaling analysis for sev- 
eral other values of V and obtained the phase diagram 
shown in the inset of Fig.l. There is a qualitative differ- 
ence between the DQMC results and the transition line 
mean field theory (MFT) predicts. Beyond V = 0.5 the 
DQMC curves rises much more sharply than the MFT 
one as a result of the competition between interplanar 
singlet formation and AFM correlation which is lacking 
in the mean-field description. Note also that known re- 
sults for the Heisenberg bilayer 15 , [l6| imply that for 
V > V c = 1.59, no order is established regardless of the 
magnitude of U. Determining the values of U c as V c is 



approached becomes quickly intractable for V > 1.0 since 
the energy scale at which magnetic correlations develop 
decreases and fluctuations in the DQMC measurements 
of long range correlations increase. 

Local moments: Within DMFT0] , the phase boundary is 
determined by a discontinuity in the double occupancy 
d, the latter being related to the local moment m by, 

m = Jj E<K) 2 > = 1 - I £<»lt»H> =1-2* 



Local moment formation is the key signature for the on- 
set of Mott insulating behavior and it has been previously 
reported to happen discontinuously at a Mott metal- 
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FIG. 2: (a) Local moment m as a function of U for 8x8 layers 
at /3 = 15 (symbols), (b) first derivative of the local moment 
with respect to U, dm/dU, shows a peak at the transition to 
the Mott insulating state, (c) A close up view of (b) shows 
the peak is no longer present for any U above V c = 1.2. 



insulator transition within other approaches as well, such 
as path-integral renormalization group [l3| and varia- 
tional Monte Carlo Figure^ shows the dependence 
of m on interaction strength for several values of V with 
no evidence of the sharp discontinuity found in DMFT. 

However, for small V (e.g. V = 0.5), we found that 
the magnetic transition is located in correspondence to 
a maximum in dm/dU (see Fig. |2Jd) . Numerical differ- 
entiation docs not allow us to establish whether there 
is an actual singularity in the behavior of dm/ dU — so 
that one could use this quantity to characterize a phase 
transition — or whether the maximum is a simple mani- 
festation of a cross over to the local moment regime. At 
larger V the value of U where the maximum appears is 
reduced (Fig. [2J; and circles in Fig.l), presumably as a 
consequence of the increased electron localization on the 
interplane bonds, whereas U c for the onset of AFLRO is 
expected to grow monotonically. This decoupling of the 
behavior of local moments from magnetism is suggestive 
of the possibility of an intervening Mott insulating state 
with no broken symmetries. 

Energy gaps and spectral functions: We now investigate 
the evolution of the energy spectra. The single-particle 
gap A sp and spin excitation gap Ag were extracted from 
the imaginary time-dependent correlations 

G w = EM T )4«r(°)>«e- rAw 

X(r) = ^(cr? CT (r)<7? CT (0))c<e-^. 

i,er 

Figure^ shows the evolution of these gaps when the in- 
teraction U is increased, for different values of the inter- 
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FIG. 3: Main panel: Gaps in the single particle spectral func- 
tion and the spin-spin correlation function are shown as func- 
tions of interaction U for several values of interlayer hopping 
V. At weak U, the spin gap As is precisely twice the single 
particle gap A sp , as expected in a Fermi liquid phase. The 
smaller panels at right show the finite size (top) and finite 
temperature (bottom) effects for the spin gap. 



layer hybridization V. Starting from the non-interacting 
limit, where the single-particle gap and the spin gap 
are expected to be respectively equal to V and 2V, the 
two quantities follow opposite evolutions regardless of 
whether there is a tendency toward AFM (small V) or 
singlet formation (large V). In particular, we found that 
for all three values of V considered in Fig. 3, the effect 
of correlation is negligible up to U ~ 2, in agreement 
with the findings of Sentef [8(. However, we observe a 
striking discrepancy between the effect of correlation in 
DMFT and in DQMC: in DQMC the single particle gap 
A sp shows no indication of the shrinking trend predicted 
by DMFT, not even in the small U limit where the role 
of long range correlation should be negligible and the 
paramagnetic solution is likely the correct ground state 
within single-site DMFT. We checked that the values of 
the gaps are converged in both size of the cluster (Fig.[3]b) 
and temperature (Fig. |3fc). 

In-layer momentum-resolved single-particle and spin 
excitation spectra (A(q, uf) and x(q,w)) are obtained by 
inverting the integral equations 



G(q,r) = / A(q 



x(q,w) 



duj. 



using the maximum entropy method [l7, 18 1. G(q, r) and 
x(q, t) are the in-layer momentum resolved counterparts 
of the correlation functions previously introduced. Fig- 
ure2]shows the single-particle (top panels) and spin (bot- 
tom panels) spectral densities and compare an AFM situ- 
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FIG. 4: Top row: single particle spectral function in the 
presence (left) and absence (right) of AFM order. Bottom row: 
same as top row but for the spin spectral function. Results 
are computed on a N = 12 x 12 cluster at ft = 20. Lines 
in the top panels are the corresponding energy bands when 
[7 = 0. 



ation (left) against ct C0/SG where no order is found by a fi- 
nite size scaling analysis (right). Thanks to the Goldstone 
theorem, the spin spectra provide a complementary indi- 
cation of the presence or lack of AFM long range order. 
We verified that, indeed, parameter regimes that were 
predicted to be AFM by scaling analysis are character- 
ized by the existence of a massless mode at (n, tt) which is 
conspicuously absent in paramagnetic cases. The single 
particle spectrum, on the other hand, helps in character- 
izing the paramagnetic phase more precisely as it shows 
flattening of the non-interacting bands. This is the typi- 
cal scenario in the presence of a Kondo resonance and it 
is therefore natural to characterize the phase at V = 1.5, 
U = 10 as a Kondo insulator. At large V, but still below 
V c , our results therefore suggest that, upon increase of U , 
the system shows a first crossover to a Kondo insulating 
state and then a transition into the anti-ferromagnet. 

Spin correlations: The real space spin correlations across 
the layers (oji ■ ayi) are shown in Fig. [5^,. The generic 
behavior of bilayer models (and related Hamiltonians like 
the periodic Anderson model) is the development of sin- 
glets with increasing V at fixed U, and the associated de- 
struction of AFLRO, signalled by a growth in (oji • 032}. 
The development of such interplane spin correlations can 
be seen in Fig. [5] by comparing the different curves at 
fixed U. The evolution at fixed V also provides consis- 
tent indications of the underlying physics previously in- 
ferred from the structure factor S* af , the local moment m, 
and the excitation gaps. Specifically, the interlayer spin 
correlations first increase as interactions are turned on, 
but then have a kink, or even turn over, as the AF phase 
is entered. For V = 0.5 for example, the kink appears at 
U « 4.0. 



The intra-plane real space nearest-neighbor spin cor- 
relations are shown in panel (b) of Fig. [5l They increase 
monotonically (in absolute value) with U for all V, indi- 
cating that the on-site Hubbard U enhances short range 
intraplane antiferromagnetism[l||. This quantity offers 
yet another local diagnostic for the onset of order as it 
shows an inflection point in close correspondence to the 
transition. Moreover, comparison of the results in the 
two panels at V = 0.5 and small U reveals that it is the 
interplane spin correlation that grows more rapidly. We 
take this as a further indication that the discrepancy in 
the behavior of the single particle gap between our cal- 
culation and single site DMFT is not due to the presence 
of intralayer short-ranged magnetic order but to the in- 
clusion, by the DQMC approach used here, of interlayer 
singlet correlations. 
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FIG. 5: Near-neighbor real space spin-spin correlations (a) 
across the layers and (b) within an individual plane. At small 
V < 1.2 the correlations converge to finite values characteriz- 
ing the magnetic ordered phase, while for large V = 2.0, the 
system is made of almost decorrelated singlets. 

Conclusions: In this paper we studied the effect of in- 
troducing local interaction in the band insulator formed 
by a bilayer with opposite sign of the hopping integral. 
We found strikingly different physics from the ionic Hub- 
bard model owing to the fact that the system is per- 
fectly homogeneous and accompanied by a tendency to- 
ward singlet formation as the band gap increases. As 
the strength of the interaction is increased, and below a 
critical interplane hybridization, a transition to a Mott 
insulator with antiferromagnetic order ensues. This tran- 
sition was studied by examining several physical observ- 
ables such as the magnetic structure factor, the local 
moments, single-particle and spin excitations resolved in 
both energy and momentum, and spin correlations. The 
behavior of dm/dU suggests that, as V grows, the mag- 
netic transition is preceeded by a cross-over into a fea- 
tureless insulating state that we characterized as a Kondo 
insulator. A more subtle question is whether such cross- 
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over may, in fact, be a transition. Obviously, this is a 
delicate point that requires validation from calculations 
on larger clusters and the use of a direct estimator for 
dm/dU rather than the finite difference employed in this 
work. In contrast to what happens in the related ionic 
model and unlike the scenario depicted by DMFT for the 
same model presently considered, the single particle gap 
increases monotonically upon increase in U and no sign 
of intermediate metallic phase is observed. 
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